Phase Coexistence of a Stockmayer Fluid in an Applied Field 



Mark J. Stevens and Gary S. Grest 
Corporate Research Science Laboratories, Exxon Research and Engineering Company, Annandale, 

NJ 08801 
(February 8, 2008) 

Abstract 

We examine two aspects of Stockmayer fluids wliich consists of point dipoles 
that additionally interact via an attractive Lennard-Jones potential. We per- 
form Monte Carlo simulations to examine the effect of an applied field on 
^ . the liquid-gas phase coexistence and show that a magnetic fluid phase does 

■ exist in the absence of an applied field. As part of the search for the mag- 
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netic fluid phase, we perform Gibbs ensemble simulations to determine phase 
coexistence curves at large dipole moments, /x. The critical temperature is 
found to depend linearly on fp' for intermediate values of /u beyond the initial 
nonlinear behavior near /u = and less than the /i where no liquid-gas phase 
coexistence has been found. For phase coexistence in an applied field, the 
critical temperatures as a function of the applied field for two different /i are 
mapped onto a single curve. The critical densities hardly change as a func- 
tion of applied field. We also verify that in an applied field the liquid droplets 
within the two phase coexistence region become elongated in the direction of 
the field. 
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I. INTRODUCTION 



Dipolar fluids have a simple, but anisotropic, long-range pair potential which presents an 
interesting new set of problems for statistical physics. Manifestations of dipolar fluids include 
ferrofluids (FF), magnetorheological fluids (MR), electrorheological fluids (ER) and polar 
fluids Hlj-Q. These fluids have great technological promise and in some cases applications, 
as they possess the dual properties of being fluid and magnetic (electric) . In most cases, an 
applied fleld is present. For the MR and ER fluids, an applied fleld is essential, since only 
induced dipoles exists. 

Simulations of simple dipolar fluids have discovered several interesting phenomena. The 
fluid can become magnetized in the absence of an applied fleld at high densities [^,^. For the 
simplest model of hard or soft sphere dipoles in zero fleld, no liquid-gas phase coexistence 
has been found in contrast to predictions of mean fleld theory In an applied fleld phase 

coexistence does occur for this system, yet it is not the usual gas-liquid or even isotropic- 
magnetic liquid phase coexistence . Coexistence in the absence of an applied fleld does not 
occur in part because the only attractive interaction is the anisotropic dipolar interaction 
which tends to produces chains instead of droplets. Phase coexistence does occur in zero fleld 
if a sufficiently strong, short-range attractive interaction is added to the dipolar interaction 

i- 

One model that exhibits two phase coexistence in zero external fleld is the Stockmayer 
fluid (SF) This system consists of long-range dipoles that have an additional short- 



range Lennard- Jones (LJ) interaction. Since the LJ interaction alone is sufficient to produce 
a gas-liquid coexistence 0] , it is clear that the Stockmayer fluid must also have a gas-liquid 



coexistence at least for small dipole moments. Recent simulations JTOl have calculated the 



coexistence curves for several nonperturbative values of the dipolar moment, /i, in zero 
external magnetic fleld. For this reason the Stockmayer fluid is a good system to study the 
effects of an applied fleld on the gas-liquid phase coexistence. 

One intriguing aspect of dipolar fluids is the existence of a magnetic fluid phase in the 
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absence of an applied field. For soft sphere dipoles, simulations have found a magnetic fluid 
phase [^,0 at high densities. Early simulations on the Stockmayer fluid for relatively small 
fi found negative pressures at these high densities ||1T|. Consequently, later studies of the 
ordering in the fluid preferred the hard or soft sphere dipolar system [^,0 . We now know that 
the negative pressures occur because of the two phase coexistence present in the Stockmayer 
fluid. Prior to this work, there have been no simulations to examine whether the magnetic 
liquid phase exists for a Stockmayer fluid at densities above the liquid coexistence curve. 
We report here results of simulations for the Stockmayer fluid that search for the magnetic 
liquid phase. The phase coexistence curves were calculated for larger values of fi where the 
magnetic liquid phase is more likely. Using these calculated coexistence curves we examined 
the regime above the liquid coexistence densities and find the magnetic liquid phase. 

The Stockmayer fluid provides a simple model system for ferrofluids Recent work 
|8|j9|JT^ has shown that a dipole with a purely repulsive core potential does not model some 
aspects of ferrofluids and that the presence of an additional attractive interaction is essential 
for phase coexistence. Experiments on hydrocarbon based ferrofluids have shown that phase 
coexistence depends strongly on the solvent |T^JT^. For some solvents phase coexistence in 
the absence of an applied field does occur []T^, while for others it does not |T^. Given that 
the same behavior occurs in simple dipolar systems as a function of the strength of the 
attractive part of the interaction, an examination of such systems is warranted. 

In the next section, we describe the simulation methods employed. We also show that 
the variant of the Stockmayer potential used in Ref. P| can be mapped to the standard 
Stockmayer potential. In section [IT| , we present calculations of the coexistence curves for 
large fi and verify the mapping determined in section |I[ The existence of the magnetic fluid 
phase is demonstrated in section 0. In section |V| we examine the effect of an applied field 
on phase coexistence. 
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II. SIMULATION METHODS 



The dipole-dipole interaction between particles i and j is 



(1) 



where /Xj is the dipole moment of the ith particle, is the displacement vector between 
the two particles and a hat signifies a unit vector. There are three relevant parameters 
for the dipolar systems: the density, p, the dimensionless dipolar coupling strength given 
by the ratio, A = fi'^/a^T, where a is the particle diameter, and the dimensionless applied 
field, rj = fiH/T. The temperature may also be an independent variable depending on the 
nondipolar interactions. In any case, the dimensionless temperature is r = 1/A. 

In addition to the dipole interaction the Stockmayer fiuid possesses a Lennard- Jones pair 
potential which models the van der Waals interactions in ferrofiuids 



Uuirij) 
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(2) 

The temperature and A are now independent quantities. The LJ system with no dipole 
interaction (/i = 0) has a critical point at T* = 1.316 and p* = 0.304 [|T^, where variables 
are given in reduced form: T* = T/e, p* = pa^, p*^ = p^ /ea^ and H* = H^a^/e. This 
critical point value is for a LJ cutoff equal to half the box length which is shifted slightly 
from the critical point of the LJ with no cutoff ||TB|. To maintain consistency we always use 
a cutoff equal to half the box length as was done in previous simulations for the Stockmayer 



fluid in zero fleld for several values of p*"^ < 6 [10|; at p* = 2, the critical point is T* = 2.09 



and p* = 0.289 |I0 



Recent simulations examined ostensibly a different system in which the strength of the 
attractive part of the LJ interaction was varied [Q. Speciflcally, they used the potential 



Ae 



(3) 



where ee is the constant used to vary the strength of the attractive interaction. We will 
refer to the system of dipoles with Uljq as the SF6 system. They found phase coexistence 



for ee ^ 0.30. For smaller values of eg, two phase coexistence was not detected by the 
Gibbs ensemble simulations and near the expected critical points chain formation was found 
instead. Thus, there seems to be a minimum amount of attraction necessary for a dipolar 
liquid phase to exist. This is consistent with the absence of phase coexistence for the soft 
sphere dipolar system, = [^[H . 



The significant feature of the interaction is not the absolute strength of the nondipolar 
attraction, but the relative strength in comparison with dipolar interaction. This is because 
the dipole interaction induces a condensed phase different from the LJ interaction, and 
consequently, the phase diagrams for LJ and strong dipolar (e.g. soft sphere dipoles) systems 
are different 0. The dipolar interaction because of its orientation dependence prefers to 
aggregate particles in anisotropic chain structures instead of isotropic droplets like the LJ 
interaction prefers. That the relative strengths are the essential quantity is evident from 
that fact that Ulm can be mapped onto Ulj [0. Mapping 



^' = (4) 

and 

e' = eel (5) 

converts the LJ6 parameters, a and e, into LJ parameters a' and e'. The density and 
temperature are mapped in the following manner 

P' = pI^T. (6) 

and 

T = T/el (7) 
We can map the dipolar system, SF6, to SF by 

P' = p/eT. (8) 
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Thus, reducing ee as was done in Ref. is equivalent to increasing the effective dipole 
moment. We confirm this mapping in section From the results of ref. p|, this mapping 
implies that for the Stockmayer fluid there is no phase coexistence for /i* ^ 5. 

Our simulations methods follow that of previous works on dipolar systems P,P, p^ , |20| , |21 



We performed simulations in zero fleld at /i* = 2.5, 3.0, 3.5 and 4.0. Simulations as a 
function of the applied fleld were performed at /i* = 1 and 2.5. The Ewald sum is used 
to evaluate the dipole interaction in periodic boundary conditions with the convergence 
parameter, a = 5.75, and reciprocal vectors are summed to lOvr/L where L is the length of 
the simulation cell. The Ewald sum contains a boundary term for the total potential energy 

where M = /^i is the magnetization. We use fiBC = making the above term zero. This 
allows the uniformly magnetized state to occur. A value of /i^c = 1 will lead to magnetic 
domains because this condition prevents the system from having a net magnetization 0. 
This form of the Ewald sum treats a periodic lattice enclosed in a large sphere or cube 
P^ . These boundary conditions are consistent with experiments performed between parallel 
plates. As mentioned above, the LJ cutoff was set at half the box length and the usual long 
range corrections were employed p3| . 



In recent density functional calculations ||2J|, it was found that the phase coexistence 



depends on shape of the system volume. In these calculations, the system was taken as an 
ellipsoid of revolution with axis of revolution k times the other two axes. The magnetic 
fluid phase is favored for large k. Simulations with the Ewald sum usually treat a cubic. 



or equivalently spherical, system [Q. The cube has sidelength, S ^ L, where L is the 
simulation box length. The simulation is generally viewed as modeling a microscopic piece 
of the bulk system. For example, with respect to ferrofluids which are often experimentally 
studied conflned between parallel plates, the simulation models a cube within the middle of 
the system. Ideally, one wants to choose fisc equal to the permeability of the bulk system, 
fJ'buik, but that requires a priori knowledge which one does not have. Recent simulations 
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have shown that the difference between using fiBC = oo and Hbc = f^buik is small as long as 
fibuik is large [g. 

To examine the liquid structure we performed constant volume canonical ensemble sim- 
ulations with N = 256 particles. These simulations ran for a least 10^ MC cycles with each 
cycle comprised of an attempt to translate and rotate each particle. The particle step size 
was chosen to achieve about 50% acceptance. The rotation step was usually set so that the 
component of the new moment in the direction of the old moment was at most 0.2. This was 
done for most cases since the acceptance rate could not be lowered to 50%. The rotation 
step size was lowered below 0.2, when that would yield a 50% acceptance. 

Gibbs ensemble simulations were performed to determine the coexistence curves and 



critical points p5| , |26[| . Below the critical temperature the Gibbs ensemble directly gives 
the two coexistence densities at a given temperature. For most of the Gibbs simulations, 
= 512. Some simulations were performed with N = 1024 to test the dependence. 
The Gibbs simulations required at least 10'^ cycles. The larger fi* simulations required runs 
2-5 times longer. Long runs were also performed for temperatures near T^. Each cycle 
included an attempt to move each particle once, 100 attempts to change the cell volume 
and 500 attempts to exchange particles between the two cells. Only a small percentage 
of exchanges is accepted in these Gibbs simulations. Thus, an extremely large number 
of attempts is needed to gather reasonable statistics. For our simulations at least 10,000 
accepted exchanges occured in 10,000 cycles. We have verified that both the pressures and 
the chemical potentials are the same in the two cells. The pressure is calculated from the 
virial expression and the chemical potential is calculated from the overlap of the particle 
insertion and extraction energy distributions . 



Our simulations in the Gibbs ensemble were performed with H* ranging from to 5. 
The critical temperature, Tc, and density, pc, were determined by fitting the calculated 
coexistence curves using the law of rectilinear diameters and the usual scaling law for the 
density with exponent (3 = 0.32 



We characterize the system structure through the order parameter. Pi, defined as 



P^= y jj, d = —M d, (10) 

where d is the director and M is the total magnetization of the system. For a completely 
magnetized system, Pi = 1. 



III. PHASE COEXISTENCE AT LARGE fi* 

In order to examine the possibility of a magnetic fluid phase, we performed simulations 
at larger fi* than previous works [|l^ . Using Gibbs ensemble simulations we determined the 
coexistence curves for /i* = 2.5, 3.0, 3.5 and 4.0. The critical points are given in Tables 
|I[ We then performed constant volume canonical ensemble simulations at densities above 
the liquid coexistence densities looking for a magnetic fluid phase. Before discussing the 
magnetic fluid phase, we discuss the fi* dependence of the critical point. 

We plot T* as a function of /x*^ in Fig. [I|(a). The open squares are from previous 
simulations []TU[, our data are the open circles and the solid circles are the mapped data of 



Ref. [§]. For /x*^ > 2 a linear fit to T* with slope 0.254 and intercept 1.06 fits the data well. 
This fit must end near /x*^ ~ 24 where no coexistence was found [§. The parameterization 
given by van Leeuwen [|l3](c) fits the T data except at large fi*"^. Clearly, the data is 
consistent with the mapping given in Eqs. For small /i*, the dipole interaction is 

effectively a van der Waals attraction [|1^. This increases the total attraction and 
consequently, Tc is raised and pc is lowered as seen in Fig. |l|. For large /i*, this basic trend 
continues and has a simple dependence. 

Similar results are found for pc in Fig. |l](b). There is more uncertainty in this data, but 
a least squares fit calculated as above gives the general decreasing trend with increasing p*. 
Here, the parameterization of van Leeuwen only works for the original data fitted. 

As noted earlier, we can also define the dimensionless (dipolar) temperature, r. The 
critical value, Tc, is more relevant to determining the dipolar structure. In Fig. |^ we plot 
as a function of p*"^. The solid line is the equation 
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r. = 0.254+1^ (11) 

derived from the linear fit in Fig. 0. The equation suggests an apparent saturation of Tc to 
about 0.25. However, the largest value of /i*^ for which phase coexistence has been found 
= 24.34 from eg = 0.30 [|) has = 0.30. Beyond /i*^ ~ 24 Eq. |ll] does not hold and 
the lowest critical dimensionless temperature is = 0.30. At higher /i* the fluid structure in 



the vicinity of the Tc from Eq. |IT] exhibits chain formation [§] . This change in fluid structure 
precludes liquid-gas phase coexistence. 



IV. MAGNETIC FLUID PHASE 

We now address the question of the existence of a magnetic fluid phase in the absence 
of an applied field in the Stockmayer fluid. A magnetic fluid phase has already been found 
in the soft sphere dipolar system for A ^ 4 Pj6|j9|,p!^ at high densities. To determine if 
the magnetic fluid phase exists for the Stockmayer fluid, we examine densities larger than 
the liquid coexistence densities, pi. We use the coexistence curves calculated in the last 
section to determine the density region of interest. The coexistence curves for larger /i* were 
calculated in the last section because a sufficiently large A ^ 1 is necessary for existence of 
the magnetic fluid phase as found for the soft sphere system. For A ^ 1, thermal interactions 
dominate and there is no magnetic fluid phase. 

We find that for the Stockmayer fiuid the structure of liquid coexisting phase in our Gibbs 
simulations is isotropic. It is possible for the coexisting liquid phase to be magnetically 
ordered p^ , p9| -pl[], but that is more likely to occur in the case of a liquid-hquid coexistence, 
where one liquid phase is isotropic and the other is magnetic. In the canonical ensemble 
simulations at T << T^, we do find some ordering, but this is most likely a finite size 
effect. Figure |^ shows the Gibbs ensemble data and the fit to the coexistence curve for 
/i* = 2.5. Obviously, far from T*, the coexistence curve will probably not follow the simple 
fit function used, but we are only concerned about having a guide to choose where to perform 
simulations. The squares show the p* and T* at which simulations were performed in this 

9 



case. A solid square denotes a negative pressure and is within the coexistence region as 
expected. Similar simulations were performed for /x* = 3.0 and 3.5. The results are given in 
Table |m|. 

At Tc, the values of A = Ac are too low for the magnetic fluid phase based on our 
experience with the soft sphere dipolar fluid. For soft sphere dipoles, we found that for 
A = 4, the magnetic fluid phase occurs close to p* = 1.0 |l^. For the Stockmayer fluid, we 



find for /i* = 2.5, Ac = 2.38. Thus, any magnetic fluid behavior, if it occurs at all, will occur 
at T much lower than Tc. The situation does not improve much with increasing fi* since the 
value of Tc = 1/Ac (see Fig. ||) saturates for large fi*. Thus, even for /x* = 4.0, Ac is only 3.17. 
To have a chance of finding a magnetic fluid phase, T must be much lower than Tc and of 
course, p > p£. As T decreases, pe increases and we will encounter the solid phase certainly 
by p* = \/2, the close packed density. For hard sphere dipoles, the liquid-solid transition for 
/i* = 2.5 occurs at p* ~ 1.0 [Q. Thus, the existence of a magnetic fluid phase depends on 
part on where the triple point temperature and density are. If the triple point temperature 
is sufficiently low, a magnetic fluid phase may exist in the region above the triple point. 
Calculated phase diagrams P3,|5D| have shown the critical point for the isotropic-magnetic 



coexistence to be at a higher temperature than the gas-liquid Tc even for p* = 2. This is 



inconsistent with our simulation results, although the density functional calculations ||24 
are primarily for ellipsiodal geometries, not the spherical geometries used in the simulations. 
The anisotropy of the ellipsoid promotes the magnetic fluid phase. For spherical geometry, 
the results are consistent with our results. Our simulations also show that the value of A 



near Tc is too small for ordering to occur. Table |ITl| shows that Pi is small for p* = 2.5 
at p* = 0.8 and T* = 2.5 which is slightly below T*. Furthermore, as T* decreases at this 
density, it is not until T* = 1.5 that Pi > 0.5, and at T* = 2.0 Pi < 1/2 for p* < 1.0. 
Thus, for the spherical geometry, the simulations show that the gas-liquid Tc is too high for 
a magnetic fluid phase to exist above it. 

For p* = 2.5, Table [IT^ shows that for T* < 1.5, Pi > 0.5. We take the magnetic hquid 
transition density to be where Pi = 0.50. There will be of course finite size effects ||3^ which 

10 



tend with increasing to lower Pi below the transition and increase it above the transition. 
Here, we are mainly interested in the existence of the transition as opposed to pinning down 
the transition point which would require simulations with much larger values of A^. For 
A = 4.19, we find the transition at p* ~ 0.90. There is the possibility that the system is a 
supercooled liquid at this density. However, since the transition for hard sphere dipoles at 
A = 6.25 is at p* ~ 1.0, we are most likely below the liquid-solid transition. Furthermore, 
we have found that the fee crystal phase melts at p* = 1.0, although the fee crystal is most 
likely not the solid phase ground state In general. Table |T| shows that for A ^ 4.0 and 
p* ^ 0.90, Pi > 0.50 and a magnetic liquid regime exists. 



V. PHASE COEXISTENCE IN AN APPLIED FIELD 

Phase coexistence in an applied field was studied at two dipole moments, p* = 1.0 and 
2.5. We chose p* = 1 since at Tc{H = 0), A = 0.71 so that the thermal and LJ interactions 
are about of equal strength. For p* = 2.5, A = 2.38 at Tc{H = 0) and the dipole interactions 
are significantly larger than the LJ interactions. This is especially true near Tc. One might 
expect that the results for the two dipoles strongly differ, but we find that much of the results 
can be described in terms of dimensionless quantities that remove the p* dependence. 

For p* = 1 and large H*, we encountered some difficulties in obtaining accurate values 
of the coexisting densities close to the critical point. This is a problem that has been 
observed before ||25|. The free energy surface becomes rather flat and the simulation can 
become trapped away from the two minima. One might expect the problem to be worse for 
p* = 2.5, but we found this not to be the case. 

Figure H shows the phase coexistence curves for p* = 1 and p* = 2.5 at selected fields. 
As H increases, the dipoles become progessively aligned with the field direction. The dipo- 
lar interaction between pairs of particles is then stronger and the dipole moments are more 
correlated. The critical temperature increases with field due to the stronger dipolar inter- 
actions as is found in the absence of an applied field. In contrast to the varying T^, the 
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critical density changes at most only slightly. For /i* = 1.0 the shape of the coexistence 
curve changes slightly as can be seen by the fact that the midpoint line becomes almost 
vertical. However, we can still fit the curves with (3 = 0.32 and the /i* = 2.5 midpoint lines 
maintain a negative slope. Because the midpoint line was more vertical for /i* = 1 than 
fi* = 2.5, more data was needed near which tended to encounter the convergence problem 
mentioned above. 

The field dependence of the critical temperature can be simplified by examining the 
critical temperature at H, Tc{H), relative to the zero field critical temperature, Tc(0) in 
terms of the ratio 

- T,(0) • ^''^ 

A collapse of the data (Fig. |^) for our two dipole moments is obtained when we plot 
versus the dimensionless field using as the temperature, 

VH = (13) 

We have drawn a least squares fit to the data excluding the H = point which fits the data 
within the uncertainty. Near H = there must be some nonlinear behavior. At large H, i.e. 
beyond the saturation field, a fi dependence is expected because saturation is /x dependent. 
We have not reached the saturation H in our simulations. In order to determine Tc for the 
infinite field, we performed Gibbs ensemble simulations with the dipole moments fixed in 
the 2;-direction. For this case, we find T* = 3.64 for /i* = 2.5. This yield = 1.38 which 
from the fit occurs for rjH = 20.5 which is beyond the range we have studied. Thus, for the 
range of H we studied, we see no effects of field saturation. 

We can also get an idea of how near the systems are to saturation by calculating the 
magnetization, M, or equivalently. Pi. In Tables [V| and |V] we list the values of the order 
parameter. Pi. We want to calculate M at the critical point. The magnetization in the 
gas phase is basically constant as the dipoles only interact weakly. The interaction energy 
between a dipole pair will be fc^T at a separation of r = A^/^cr. For p* ~ 0.30, the average 
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separation is a = p~^/^ ^ 1.5cr. We are working in the range of A < 4, which gives r = 1.6cr. 
Thus, the dipolar interaction energies are at most equal to ksT. We take = M(Tc) to 
be the value of M for the gas phase. In Fig. ^ we plot M/ fiN for the two cases. Neither 
case has reached saturation, as we expected since there is no /i dependence in our . 

One of the important question concerning dipolar fluids (FF, MR, and ER) presently 
under consideration PP,p!^,p^-p6| is the shape of the liquid coexisting phase. We can only 
examine the structure at the particle level and are mainly concerned with confirming the 
basic effect of elongation of a liquid droplet in the applied field. We have performed a 
constant volume canonical ensemble simulation for /i* = 2 at p* = 0.01, T* = 1.0 and 
H* = and 1.0. In Fig. ^we show projection plots for the two different H*. In the absence 
of a field, the droplets are spherical on average as expected, and in the presence of a field 
the droplets become extended along the field direction {z in Figure) and in this case two of 
the droplets at H* = coalesced into one. Thus, we find the two main effects of an applied 
field on liquid droplets: elongation in the field direction and coalescence . 



VI. CONCLUSION 

The results of the the present simulations give a better understanding of the phase 
diagram of the Stockmayer fluid in both zero or nonzero fields. The mapping of the SF6 
system onto the SF system shows that dipole interaction strength in comparison with the 
LJ interaction strength is the determining quantity for the occurance of phase coexistence. 
The mapping is particularly useful as discussed below when comparing with experimental 
systems as discussed below. The /i dependence of the critical point is rather simple for most 
of the range over which SF phase coexistence occurs. In this range, the temperature, T* 
depends linearly on p^. Linear behavior is also true for pc, although there is more uncertainty 
here. The magnetic fluid phase does appear at high densities for sufficently large /i. At least 
for the spherical geometry, the tricritical temperature appears to be below T^. In an applied 
fields, a p-independent scaling can be obtained at least for a broad range of fields for the 
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critical temperature as a function of applied field when the temperature is scaled by the zero 
field Tc. 

One of the important conclusions of previous works PJT^ is that the hard or soft sphere 
dipolar systems were insufficient as models for ferrofiuids (and most likely MR and ER 
fluids). Some added central force attraction is required such as found in the SF or SF6 
potentials. Dipole moments in terms of A are about 1, although ferrofiuids possess a large 
polydispersity However, the residual strength of an attractive interaction such as van 
der Waals is unknown. Thus, the value of ee must be determined from some experimental 
data. If phase coexistence occurs in zero field, one way of determing ee is to map the value 
of Tc for H = 0, and the given value of fi onto the plot of versus eg. In the cases where no 
phase coexistence occurs at H = 0, one must determine eg from the value of at H > 0. 
The linear relationships between or pc and /i^ simplify this procedure somewhat. In some 
of the experiments [0] phase coexistence curves are not fully measured, but slices in the 
H-p plane at fixed T are measured. To compare with these experiments requires calculation 
of coexistence curves not only at several H, but also several T and eg in order to determine 
the correct eg and the make the correct slices. 

The structure at the particle level in ferrofiuids has not been resolved experimentally. 
Simulations naturally offer a means to examine this structure. The structure of the columnar 
objects |jl|j3^ that form in ferrofiuids is an open question. In the simulations of soft sphere 
dipoles, chains form consisting of connected single particles as in a polymer P,p^. These 
chains do not aggregate to form the columnar structures as has been observed experimentally. 
In contrast, we find in the simulations of the Stockmayer fluid that the polymeric chains 
do not form. In an applied field, the zero field droplets distort becoming elongated in the 
field direction. Droplets are also seen to coalesce. This suggests that in ferrofiuids that the 
columns are formed by initial formation of liquid droplets that are distorted into an elliptical 
shape and coalesce in an applied field to form columns spanning the experimental cell. 
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TABLES 



TABLE L Critical points. 


* 


H* 


C 


Pc 


2.5 


0.0 


2.63(1) 


0.29(1) 


3.0 


0.0 


3.35(1) 


0.25(1) 


3.5 


0.0 


4.20(1) 


0.24(1) 


4.0 


0.0 


5.07(5) 


0.24(1) 


LO^ 


0.0 


1.41(1) 


0.30(1) 


1.0 


1.0 


1.44(1) 


0.32(1) 


1.0 


2.0 


1.49(1) 


0.33(1) 


1.0 


3.0 


1.51(1) 


0.32(1) 


2.5 


0.5 


2.71(1) 


0.285(1) 


2.5 


1.0 


2.78(1) 


0.285(1) 


2.5 


2.0 


2.89(1) 


0.302(1) 


2.5 


5.0 


3.15(1) 


0.278(1) 


2.5 


oo 


3.64(1) 


0.303(1) 



■Data from Ref. [|. 
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TABLE II. Phase coexistence data for fx* = 3.0, 3.5 and 4.0 









p* 

9 




T3* 


2.70 


3.0 


0.011(1) 


0.021(1) 


0.68(1) 


0.01(5) 


2.85 


3.0 


0.020(2) 


0.036(3) 


0.63(2) 


0.02(5) 


2.95 


3.0 


0.033(5) 


0.053(3) 


0.59(2) 


0.05(2) 


3.00 


3.0 


0.060(9) 


0.057(5) 


0.57(2) 


0.05(6) 


3.10 


3.0 


0.060(8) 


0.077(1) 


0.57(2) 


0.09(6) 


3.20 


3.0 


0.072(9) 


0.099(5) 


0.46(3) 


0.09(1) 


3.25 


3.0 


0.086(14) 


0.099(2) 


0.43(4) 


0.11(1) 


3.30 


3.0 


0.13(1) 


0.124(4) 


0.39(2) 


0.12(1) 


3.35 


3.0 


0.21(2) 


0.148(5) 


0.30(2) 


0.14(2) 


3.50 


3.5 


0.011(1) 


0.025(2) 


0.64(2) 


0.03(9) 


3.70 


3.5 


0.017(1) 


0.043(2) 


0.58(2) 


0.03(3) 


3.80 


3.5 


0.031(8) 


0.046(6) 


0.54(2) 


0.03(2) 


3.90 


3.5 


0.053(6) 


0.070(6) 


0.50(4) 


0.06(3) 


4.00 


3.5 


0.06(1) 


0.084(4) 


0.46(2) 


0.09(3) 


4.05 


3.5 


0.065(5) 


0.091(9) 


0.45(2) 


0.10(2) 


4.10 


3.5 


0.08(2) 


0.101(4) 


0.41(4) 


0.12(2) 


1.15 


3.5 


0.10(2) 


0.111(1) 


0.39(1) 


0.11(3) 


4.70 


4.0 


0.03(6) 


0.04(5) 


0.48(8) 


0.01(4) 


4.80 


4.0 


0.05(2) 


0.07(2) 


0.53(2) 


0.10(3) 


4.90 


4.0 


0.05(1) 


0.072(4) 


0.50(2) 


0.07(4) 


4.95 


4.0 


0.07(1) 


0.09(1) 


0.44(2) 


0.08(3) 


5.00 


4.0 


0.06(1) 


0.08(1) 


0.41(2) 


0.07(3) 


5.02 


4.0 


0.14(2) 


0.08(5) 


0.41(4) 


0.09(2) 
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TABLE III. Data for dense liquid phase for TV = 256. 





As 




P* 


p* 


Pi 


2.5 


2.50 


2.5 


0.80 


2.83(3) 


0.08(1) 


2.5 


3.12 


2.0 


0.75 


0.33(1) 


0.09(1) 


2.5 


3.12 


2.0 


0.80 


0.98(4) 


0.11(4) 


2.5 


3.12 


2.0 


0.85 


1.97(3) 


0.13(3) 


2.5 


3.12 


2.0 


0.95 


5.21(2) 


0.21(3) 


2.5 


4.16 


1.5 


0.80 


-0.93(5) 


0.24(2) 


2.5 


4.16 


1.5 


0.85 


-0.41(4) 


0.36(9) 


2.5 


4.16 


1.5 


0.90 


0.36(4) 


0.51(2) 


2.5 


4.16 


1.5 


0.95 


1.57(7) 


0.62(1) 


2.5 


4.16 


1.5 


1.00 


3.26(5) 


0.69(2) 


3.0 


4.09 


2.2 


0.80 


-0.35(4) 


0.17(10) 


3.0 


4.09 


2.2 


0.90 


1.74(7) 


0.33(5) 


3.0 


4.50 


2.0 


0.90 


0.60(15) 


0.55(6) 


3.0 


4.50 


2.00 


1.00 


3.42(9) 


0.71(1) 


3.5 


4.08 


3.0 


0.80 


0.58(5) 


0.20(4) 


3.5 


4.08 


3.00 


0.90 


2.70(7) 


0.29(5) 


3.5 


4.45 


2.75 


0.80 


0.00(14) 


0.22(8) 


3.5 


4.45 


2.75 


0.90 


1.62(12) 


0.41(10) 


3.5 


4.90 


2.50 


0.90 


0.37(17) 


0.64(2) 


3.5 


4.90 


2.50 


1.00 


2.98(7) 


0.76(2) 
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TABLE IV. Phase coexistence data in fields for /x* = 1.0 



T* 


H* 


Pi 


P* 

g 


p(9) 


Pi 


P* 


^1 


1.10 


1.0 


0.025(2) 


0.023(2) 


0.24(4) 


0.724(5) 


0.024(7) 


0.477(3) 


1.20 


1.0 


0.047(3) 


0.043(3) 


0.21(2) 


0.669(5) 


0.046(7) 


0.43(1) 


l.,30 


1.0 


0.080(5) 


0.071(3) 


0.22(5) 


0.602(5) 


0.072(6) 


0.366(7) 


1.35 


1.0 


0.113(7) 


0.092(3) 


0.21(7) 


0.554(7) 


0.083(8) 


0.341(7) 


1.37 


1.0 


0.122(9) 


0.098(4) 


0.22(2) 


0.54(1) 


0.09(2) 


0.31(1) 


1.40 


1.0 


0.15(1) 


0.11(3) 


0.22(2) 


0.48(2) 


0.100(7) 


0.30(3) 


1.30 


2.0 


0.077(6) 


0.069(2) 


0.46(1) 


0.629(6) 


0.067(7) 


0.60(1) 


1.35 


2.0 


0.090(5) 


0.079(4) 


0.45(2) 


0.590(4) 


0.08(2) 


0.57(1) 


1.40 


2.0 


0.109(5) 


0.096(2) 


0.44(1) 


0.544(8) 


0.094(9) 


0.55(1) 


1.44 


2.0 


0.157(9) 


0.117(2) 


0.45(2) 


0.50(1) 


0.118(6) 


0.53(2) 


1.46 


2.0 


0.186(8) 


0.129(3) 


0.45(1) 


0.49(2) 


0.132(7) 


0.52(1) 


1.48 


2.0 


0.21(3) 


0.138(5) 


0.46(1) 


0.45(2) 


0.141(8) 


0.51(1) 


1.30 


3.0 


0.056(8j 


0.055(5) 


0.60(1) 


O.D46(3j 


0.05(1) 


0.70(1) 


1.35 


3.0 


0.078(4) 


0.072(2) 


0.60(1) 


0.614(8) 


0.06(1) 


0.69(1) 


1.40 


3.0 


0.099(5) 


0.088(2) 


0.58(1) 


0.575(17) 


0.08(2) 


0.67(1) 


1.45 


3.0 


0.132(11) 


0.108(6) 


0.60(1) 


0.529(6) 


0.104(6) 


0.65(1) 


1.45 


3.0 


0.129(11) 


0.108(6) 


0.58(1) 


0.518(16) 


0.103(9) 


0.65(1) 


1.47 


3.0 


0.157(6) 


0.120(1) 


0.58(1) 


0.51(2) 


0.12(2) 


0.65(1) 
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TABLE V. Phase coexistence data in fields for /x* = 2.5 



T* 


H* 


Pi 


P* 
a 


p(s) 


Pi 






2.30 


0.5 


0.03(1) 


0.045(4) 


0.20(1) 


0.66(1) 


0.05(1) 


0.53(2) 


2.40 


0.5 


0.04(1) 


0.06(1) 


0.16(1) 


0.62(1) 


0.07(2) 


0.50(2) 


2.50 


0.5 


0.09(1) 


0.092(3) 


0.21(2) 


0.56(1) 


0.09(3) 


0.44(3) 


2.60 


0.5 


0.16(3) 


0.118(6) 


0.23(2) 


0.49(2) 


0.12(1) 


0.40(1) 


2.65 


0.5 


0.14(1) 


0.14(1) 


0.20(5) 


0.44(4) 


0.14(2) 


0.36(1) 


2.70 


0.5 


0.21(5) 


0.16(1) 


0.22(6) 


0.38(5) 


0.16(3) 


0.33(1) 


2.40 


1.0 


0.030(5) 


0.048(6) 


0.34(2) 


0.66(1) 


0.05(3) 


0.68(1) 


2.50 


1.0 


0.046(5) 


0.069(4) 


0.33(4) 


0.61(1) 


0.07(2) 


0.64(1) 


2.65 


1.0 


0.09(2) 


0.109(4) 


0.38(4) 


0.51(2) 


0.10(4) 


0.59(1) 


2.70 


1.0 


0.11(2) 


0.127(8) 


0.39(1) 


0.49(3) 


0.136(3) 


0.57(1) 


2.72 


1.0 


0.13(2) 


0.133(7) 


0.40(2) 


0.47(1) 


0.13(1) 


0.56(2) 


2.73 


1.0 


0.12(2) 


0.133(10) 


0.39(4) 


0.44(3) 


0.13(1) 


0.54(2) 


2.74 


1.0 


0.15(2) 


0.137(7) 


0.43(4) 


0.45(1) 


0.14(1) 


0.55(1) 


2.40 


2.0 


0.019(1) 


0.034(4) 


0.61(4) 


0.702(3) 


0.03(5) 


0.80(1) 


2.50 


2.0 


0.028(3) 


0.047(3) 


0.59(3) 


0.664(5) 


0.06(3) 


0.78(1) 


2.60 


2.0 


0.040(6) 


0.066(6) 


0.59(3) 


0.623(6) 


0.07(6) 


0.77(1) 


2.70 


2.0 


0.069(11) 


0.093(6) 


0.60(2) 


0.585(13) 


0.10(1) 


0.75(1) 


2.80 


2.0 


0.10(2) 


0.118(10) 


0.60(2) 


0.516(6) 


0.13(1) 


0.72(1) 


2.85 


2.0 


0.13(2) 


0.136(3) 


0.61(2) 


0.50(2) 


0.14(1) 


0.72(1) 


2.80 


5.0 


O.U43(8J 


O.U7U(5) 


0.82(1) 


O.D09(Dj 


U.0d3(5) 


0.86(1) 


2.90 


5.0 


0.058(6) 


0.091(6) 


0.80(1) 


0.579(8) 


0.09(2) 


0.86(1) 


3.00 


5.0 


0.10(2) 


0.118(7) 


0.80(1) 


0.54(2) 


0.12(2) 


0.85(1) 


3.05 


5.0 


0.09(1) 


0.13(1) 


0.80(1) 


0.48(1) 


0.122(8) 


0.86(1) 


3.08 


5.0 


0.12(2) 


0.14(1) 


0.80(1) 


0.46(3) 


0.14(3) 


0.83(1) 


3.10 


5.0 


0.12(2) 


0.14(2) 


0.80(1) 


0.44(2) 


0.14(1) 


0.83(1) 
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FIGURES 

FIG. 1. In (a) the critical temperature, T*, as a function of the dipole moment squared, fi*'^, 
and in (b) the critical density p* is plotted. The open points are from previous works and the 
solid points are from this work. The solid lines are a least squares fit to the nonzero p* data. The 
dotted lines are the parameterizations given in ref. [jl^(c). The uncertainty in T* is about ±0.01 
for all data which is smaller than the points. 

FIG. 2. The dimensionless temperature, Tc, is plotted versus ^u*^ with the same point types as 
in Fig. The solid line is obtained from the least squares fit in Fig. 

FIG. 3. The coexistence curve calculated in the Gibbs ensemble simulation in the absence of an 
applied field at /i* = 2.5 is plotted. The open circles are the coexisting density points found in the 
simulations. The solid circle represents the critical point calculated along with the fitting curve as 
described in the text. The square represents points where canonical simulations were performed 
to determine the existence of the magnetic fluid phase (see also Table III). 



FIG. 4. The coexistence curves within an applied field are plotted (a) for ^u* = 1 at H* = 0, 1, 
2 and 3, and (b) for ^* = 2.5 at H* = 0, 0.5, 1.0, 2.0 and 5.0. 

FIG. 5. The critical temperature for systems exhibits some scaling when plotted versus the 
dimensionless field. In order for the data of the two p* to be on the same curve, we scale T*{H) 
by the zero field T*. This is also done in the dimensionless field. The solid line is a least squares 
fit to the nonzero field data. 

FIG. 6. The magnetization per particle versus the applied field for p* = \ (triangles) and 
/i* = 2.5 (squares). 

FIG. 7. Projection plots for (a) i7* = and (b) H* = 1.0 for p* = 2 at p* = 0.1 and T* = 1 
show the effect of an applied field on the liquid droplet shape in the coexistence region. The field 
is parallel to the z-direction. 
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